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We used various analytical and numerical techniques to elucidate signal propagation in a small enzymatic 
cascade which is subjected to external and internal noise. The nonlinear character of catalytic reactions, which 
underlie protein signal transduction cascades, renders stochastic signaling dynamics in cytosol biochemical 
networks distinct from the usual description of stochastic dynamics in gene regulatory networks. For a sim- 
• pie 2-step enzymatic cascade which underlies many important protein signaling pathways, we demonstrated 

that the commonly used techniques such as the linear noise approximation and the Langevin equation become 
inadequate when the number of proteins becomes too low. Consequently, we developed a new analytical ap- 
■ proximation, based on mixing the generating function and distribution function approaches, to the solution of 

the master equation that describes nonlinear chemical signaling kinetics for this important class of biochemical 
reactions. Our techniques work in a much wider range of protein number fluctuations than the methods used 
previously. We found that under certain conditions the burst-phase noise may be injected into the downstream 
0^ ' signaling network dynamics, resulting possibly in unusually large macroscopic fluctuations. In addition to com- 

puting first and second moments, which is the goal of commonly used analytical techniques, our new approach 
provides the full time-dependent probability distributions of the colored non-Gaussian processes in a nonlinear 
signal transduction cascade. 

Keywords: Stochastic Processes, Nonlinear Chemical Kinetics, Signal Transduction, Signal Amplification, Strong Fluctua- 
tions, Master Equation 

o 

C7" 1 I. INTRODUCTION 



y—i i Biochemical signaling networks mediate information flow in cells, regulating important cellular processes such as cell 
metabolism, motility, gene expression and cell cycle 1 ' 2 . A quantitative understanding of signal transduction is essential for 
realizing the larger goal of modeling the cell. Until recently, biochemical reaction networks were analyzed mainly with the help 
of chemical kinetics equations, exploring such issues as robustness, sensitivity, and bistability^JJiii. This deterministic descrip- 
tion, though valid in the asymptotic limit of a macroscopic system (as for a reaction in organic chemist's test-tube), fails when 
the number of reacting proteins is too small. This is often the case in the cell, which is a mesoscopic system. Therefore, instead 

\Q | of a smooth deterministic course, the fundamentally random nature of chemical reactions results in "noisy" reaction trajectories 
in individual cells. The resulting heterogeneous response of an ensemble of cells to a particular external signal 8 necessitates 
going beyond chemical kinetics, relying instead on the theory of stochastic processes&iSiAA*^ to describe signaling dynamics. 

• <-h . Previous theoretical and experimental results have strongly suggested that stochasticity is an important component in the 
' dynamical processes in a ce iiI2iI2iI2iI£iASiIiI£iA(i. For instance, in signal transduction, fluctuations induce not only quantitative 
q-< changes of threshold values 20 but also qualitative changes such as oscillations, transitions between different stable states, and 
stochastic resonances that may increase sensitivity or stability of the cellular signal processing^iiiSiSiSi 2 ^ 2 ^. The resulting large 

. £h ! variability in cell-to-cell response to the same external stimulus 8 ' 27 allows a cell colony to adapt to varying environment. 

The connection between deterministic and stochastic chemical kinetics is analogous to that between classical and quantum 
mechanics^ 2 ** In particular, the number of degrees of freedom in stochastic description of chemical kinetics is immensely 
larger compared with the deterministic description. Consequently, the equations of stochastic chemical kinetics are difficult to 
solve, and mainly numerical results have been discussed in prior work. Small sizes and intrinsic complexities of cells allow 
for the emergence of large fluctuations that are coupled with other rich dynamical processes. Thus, the resulting formidable 
complexity of these dynamical systems has hindered obtaining a complete solution to the problem of stochastic dynamics in 
biochemical reaction networks. How to construct approximate analytical solutions to understand qualitatively the complex 
signaling dynamics is a challenging problem which stands in the research frontier of non-equilibrium statistical mechanics and 
dynamical systems 20 * 3 ^. 

A number of tools have been developed to deal with the randomness in chemical reactionsi^iiSiSi^i. The Gillespie simulation 
algorithm 35 - 36 ' 37 , the Langevin equation and the Fokker-Planck equation2L2i22i22i2!iiI are among the most commonly used. These 
methods have been applied fruitfully to study signal transduction processes^iSiiiil 2 *! 3 .. However, a number of constraints limit the 
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applicability of these methods. Gillespie simulation is essentially a Monte Carlo algorithm which simulates chemical reactions 
as a series of transitions between reaction states with transition rates determined from microscopic kinetics. One simulation 
corresponds to one possible reaction trajectory. When the system is large, one or two reaction paths provide qualitative or 
even quantitative information of the system dynamics. However, to get good statistics, a large number of paths are required 
when the system size is moderate or small. This, in many cases, may be computationally expensive. Moreover, although the 
computational cost is dominated by the fastest reaction in a cascade, the time scale of interest is likely given by the slowest 
reaction, perhaps many orders of magnitude slower than the fast reaction. In this situation, commonly encountered in biological 
signaling networks, producing a single trajectory is computationally forbidding. Under special conditions these simulations can 
be accelerated 9 ' 36 - 4 '-Ili^i^ilM^, but for a general reaction network this is still an active area of research. On the other hand, 
a qualitative understanding of the reaction network dynamics is essential for learning control of biochemical processes in the 
ce jj7 j 50 jj j s difficult to extract such understanding from Gillespie simulations alone, without the guidance from a complementary 
analytical model. 

A number of approximate analytical techniques are used to solve the chemical master equation. Among the most commonly 
used ones, the linear noise approximation 38 ' 39 , is an effective weak-noise expansion based on the fact that fluctuations are of order 
VQ for a system of size f2. It works well for large systems where fluctuations are small such that the probability distribution 
is centered narrowly around deterministic orbits determined from chemical kinetics?-!.. However, molecular discreteness and 
large fluctuations in cellular biochemical reactions, combined with nonlinear effects, may generate strong correlations along a 
pathway, leading to the formation of characteristic patterns both in time and spaceii&ISiSIi&SSi. 

To account for such patterns, we developed an approach which incorporates large fluctuations, going beyond the commonly- 
used small noise and continuous approximations. We focus our efforts on a specific 2-step signaling cascade, consisting of a 
unary reaction of receptor activation upstream and a subsequent binary reaction of enzyme activation downstream. Biochemical 
signaling networks are comprised of only few signaling elements, and the 2-step signaling cascade considered in this work is 
among most fundamental building blocks. The nonlinearity of the catalytic reaction in the second step is the main source of diffi- 
culties in obtaining a comprehensive analytical solution to stochastic signal dynamics in this amplification cascade. In this work, 
we developed high-quality approximate solutions to the stochastic signal propagation dynamics in a 2-step cascade, uncovering 
that fluctuations may broaden the average chemical kinetics signals such that transient, highly non-Gaussian distributions may 
emerge, due to interplay between discrete noise and nonlinearity. We could not reproduce this effect using the chemical Langevin 
equation, since the latter is based on the continuous approximation, ignoring molecular discreteness. 

The paper is organized as follows. After commenting on the strengths and weaknesses of the traditional modeling techniques, 
we solve the master equation for the 2-step cascade with the help of generating functions. Then the new formalism is used to 
elucidate how noisy signal may be controlled in an unbranched signaling pathway. In particular, the upstream fluctuations may 
propagate downstream without dissipation, resulting in large downstream fluctuations even in the limit of a macroscopic size 
downstream signaling network. Although we focused in this work on an important, yet specific enzymatic cascade, the hybrid 
generating function - distribution function technique introduced in this work (Section HVl may be generalized to treat more 
complex biochemical pathways. To demonstrate the usefulness of this hybrid "smooth distribution" method, we consider in 
Section lTVb a self-dimerization of the receptors in the first step of the cascade activation, which enhances the nonlinear character 
of the 2-step cascade. In a forthcoming publication we will illustrate the use of time-dependent basis functions, developed in 
this work, in the variational solution of stochastic chemical kinetics equations in signaling cascades comprised of several steps 
and containing feedback loops. 



II. MASTER EQUATION TREATMENT OF REACTION TRAJECTORY REALIZATIONS 

Signal transduction often starts at the cell membrane, where external ligands, such as hormones or small peptides, bind and 
activate cell surface receptors. In turn, the activated receptors send the signal downstream, usually by phosphorylating specific 
proteins inside the cell 1 . These proteins then activate other cytosol proteins further downstream. The process goes on so that 
the signal propagates through the cell in a predetermined way. However, the signaling dynamics is not uniform when protein 
numbers are small, which is often the case in the burst phase of the cascade activation. Because of the fundamentally random 
nature of chemical reaction dynamics 31 , each trajectory realization is different from others, even when the same initial conditions 
are used. This behavior is called trajectory variability in the literature 27 . If the variation is large, then a stochastic description 
becomes necessary. 

Our long-term goal is to understand qualitatively how large signaling networks process information. The 2-step amplification 
cascade, shown in Fig. [2 can be regarded as the building block of more complex cascades. In this simple reaction scheme, 
without feedback loops, R represents an inactive receptor, which becomes activated into R* upon binding of an external ligand 
(stimulus). When the receptor is activated, it acts as an enzyme, catalyzing the phosphorylation of the next kinase downstream 
(A+R — ► A* +R*) with a rate fi. A* spontaneously decays to A with a rate A. Although the R* reaction is unary and independent 
of the A reaction, the latter one is binary, making the system nonlinear , thus, different from those usually considered in the gene 
regulatory networks^^iS. 
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FIG. 1: An inactive receptor R, when activated by a signal, activates downstream protein A. 




cGMP Hydrolysis 

closure of cGMP-gated Na + /Ca 2+ channels 



FIG. 2: A schematic depiction of the visual signal transduction cascade. 



This simple 2-step cascade is commonly embedded in the onset of a reaction pathway of many important signaling 
cascades^®. For example, the visual signal transduction pathway, which takes place in retinal rod cells, is shown in Fig. |2f 9-61 ■ 
Rhodopsin receptors (Rh), located in small discs in the outer segments of the rod cell, contain a light-sensitive molecule, reti- 
nal. Upon incidence of photons, the retinal molecules undergo isomerization, leading to a subsequent activation of rhodopsin 
receptors (Rh*). The latter act as a catalyst to replace the GDP by GTP in a G-protein, called transducin. The next enzyme in 
the cascade, the cGMP phosphodiesterase (PDE), is then activated to (PDE*) by the active G-protein GTP complex. PDE* hy- 
drolyzes the cGMP, which leads to the closure of cGMP-gated Na + /Ca 2+ ion channels and, thus, reduces the influx of Na + /Ca 2+ 
flow. The rod cell becomes hyperpolarized, resulting in less neurotransmitters being released from the synaptic terminal of the 
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rod cell. This change is immediately picked up by a second cell and passed as a signal to the central neural system. In this 
example, the initiating steps of the visual signal transduction pathway, 

Rh <-> Rh* , G-GDP + Rh* <-► G-GTP + Rh* 

are similar to our 2-step cascade. The amazing ability of retina rod cells to detect single photons, in the presence of external and 
internal noise, has been discussed in prior works 62 . Although understanding visual signal transduction is not the focus of this 
work, the formalism of large fluctuations, developed in the current work, may serve as a basis for further elucidating this issue. 

We regard the R — ► R* reaction (Fig.^ as a Poisson process with rate g, which corresponds to an abundance of R receptors 
and scarce ligand presence. For instance, the light perception in a dark roomS^ may be described as a Poissonian bombardment 
of photons on the retina cell surface. If the intensity of photons is low, there always exist a sufficient excess of inactivated photon 
receptors (R) such that the arrival events of photons dominate the process. R* spontaneously decays back to R at a rate k. The 
average number of activated receptors (R*) depends both on the incidence rate g and decay rate k. If R — > R* is considered to 
be an ordinary first order chemical reaction, instead of a Poisson process, all methods described below would still apply, with 
only minor modifications. Our current investigation is restricted to the 0-D treatment of space. This is a good approximation if 
the reaction network is confined to a small enough spatial region, having a linear dimension of £ (so-called Kuramoto length^), 
such that particles diffuse across the region fast compared to the typical reaction times. Using the reaction parameters from 
published models of the EGF/MAPK signaling cascada^^ 6 !, we estimated £ to be in the range from 0.3 \xm to 5 /im. In an 
ongoing work we are incorporating an explicit treatment of space into our analysis. 

150 




FIG. 3: Time evolution of the variance of A* computed from Eq.Q(.so/;c/ line), Langevin equation (dashed line) and Ed. 1331 (circles). For 
(g, k, n, A) = (0.2, 0.1, 0.02, 0.15) with initial condition (N R , N R * , Na, N a * ) = (100, 0, 100, 0). 

When the number of protein copies in the signaling network is small, the evolution of the averages, described by ordinary 
chemical kinetics, is inadequate to characterize the system dynamics. An example of the evolution of the A* protein number 
variance, shown in Figure [5] indicates that the amplitude of the fluctuations is about 10 in the steady state. This is significant 
when compared to the deterministic average of about 20. Clearly, a stochastic description is required in this case. The chemical 
master equation 31 is a starting point for studying large variations of individual stochastic trajectories from the average path. For 
the 2-step cascade described above (Fig.Q, we denote by P(m, n) the probability of having m copies of R* and n copies of A 
at a particular time point. The time evolution of this probability distribution is determined by the following master equation 

dP 

—j^(m,n) = fi[— mnP(m, n) + m{n + l)P(m, n + 1)] + A[— (N — n)P(m, n) 
+ (N — n+ l)P(m, n - 1)] + g[-P(m, n) + P(m - 1, n)} 

+ k[-mP(m,n) + (m + l)P(m + l,n)}, (1) 

which expresses the transition rates of probabilities from time t to time t + dt in terms of the probability distribution at time t. 
The sum of A and A* is taken to be constant throughout the reaction (N in Eq. Q). Another way to think about the chemical 
reaction dynamics given by Eq. Q is to view it as a random walk on a two-dimensional lattice of integer coordinates m and n 
with position-dependent jump probabilities^. 

Since the master equation provides a full description of stochastic chemical process, solutions of the set of coupled ordinary 
differential equations in Eq. Q provide all necessary information to analyze the signaling dynamics. However, an exact analyt- 
ical solution for the system of ODEs in Eq. Q is not known. Direct numerical integration is also difficult due to the enormous 
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number of ODEs (10 4 -10 8 ) for a cascade that may contain up to 10 2 -10 4 proteins of each species. Simulations based on the 
Gillespie algorithm may be used to estimate P(m, n) in Eq. 0, however, they are computationally inefficient and are hard to 
interpret, as discussed earlier. Thus, to gain qualitative insight into stochastic signal transduction by the 2-step cascade in Fig.[0 
it is important to develop an approximate analytical solution to Eq. Q. 

Physical considerations are crucially important in making sensible approximations to the master equation. Several numerical 
and analytical techniques are available to account for the noise effect on signaling dynamics, including the chemical Langevin 
equation 41 and van Kampen's i7-expansior. 31 . But they are most useful only when the particle numbers are large and fluctuations 
are relatively small. In Appendix A we derive an 17-expansion for the 2-step cascade Q. In the remaining part of the paper, we 
develop alternative solution schemes to solve equations like Eq. Q, based on the generating function approach, to obtain signal- 
ing dynamics in the regime of large fluctuations, where traditional analytical techniques are no longer applicable. Comparison 
will be made between our method and the results from Langevin equation or Sl-expansion. 

Generating function approaches have been used in prior works to elucidate stochastic processes in gene regulatory 
networks 55,57 . The novelty of our work is to extend these techniques to nonlinear chemical processes in protein signal am- 
plification cascades. The nonlinearity of chemical reactions greatly enriches the dynamical behavior of signaling cascades. 
However, to develop a robust analytical approach to solving stochastic amplification dynamics, substantial additional difficulties 
need to be overcome compared to describing noisy dynamics in linear biochemical networks. 



III. THE GENERATING FUNCTION APPROACH 



To treat signal transduction in a wider range of protein numbers in the cell and to analyze the effects of large fluctuations, we 
have developed a new approach, based on generating functions, to solve the master equation. A generating function encodes 
probability distributions as its Taylor series coefficients. As a result, the enormous set of ODEs in Eq. Q are compactified into 
a single PDE. Thus, the evolution of the probability distribution can be obtained by solving this one PDE for the generating 
function. Since even for medium-sized cascades, there are an astronomical number of ODEs in the master equation formalism, 
an approximate generating function greatly facilitates qualitative and quantitative analysis of strongly-fluctuating dynamics in 
biochemical reaction networks. 

As an example, for the 2-step cascade, we define a generating function through the following power series 



y 



(2) 



which satisfies the time evolution equation 

(1 - y){nx 



a* 
at 



dxdy 



d a* 

\N + \y — )y + g(x - 1)* - k(x - 1) — 
oy ox 



(3) 



Our next goal is to develop approximate techniques to solve Eq. 0. We know that the solution of Eq. Q is an analytic function 
of x , y with nonnegative time-dependent coefficients. The highest derivative in Eq. Q is d 2 ^ /dxdy which reflects the binary 
chemical reaction between R* and A. If this term is omitted, Eq. can easily be solved by the method of characteristics. 
However, this would completely alter its physical content. On the other hand, we notice that the generating function of the R* 
distribution 4>{x) = ^{x, 1) does not depend on the dynamics of A and obeys the following PDE 



which can be solved exactly, resulting in 



-= 9{x -l)-k(x-l)-, 



( a; )=exp[|(x-l)(l-e- fct )] 



(4) 



(5) 



where the initial condition Nr* = at t = was used. The R* probability distribution, P(m), is given by the coefficients of 
the series expansion of Eq. (0 



c) = ^exp(-|(i- e --))(i- e --r(fr^ 



(6) 



resulting in, 



P(m) = exp(-(l - e- kt )g/k)(l - e~ fc *) m (|) m /m! . 



(7) 
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Therefore, the time-dependent distribution of R* is Poissonian and relaxes to a stationary distribution with the rate k. Although 
this distribution is generated by both the birth and decay of R*, the relaxation is independent of birth rate g. From Eq. 0, the 
average and the variance of Nr* are also easily calculated, (Nr*) = (a 2 ) = (1 — e~ kt )g/k. 

Next, we build up the cascade by considering also the reactions involving A. In particular, we construct a series expansion^ 
in x with time-dependent functions 4> m (y), 

oo 

^(x,y)=Y J <t>m(y)x m . (8) 

m=0 

Thus, with each state containing to R*'s, we associate a distribution of A, which may be computed from cf> m (y). The new 
functions <j) m (y) satisfy 



— ^ = (1 - y)(/j.m- XN + Xy—)(j) m + km(4> m+ i - <j> m ) + g(<j> m -x - <t>m) + &<?Wi ■ (9) 

at ay ay 

Here an infinite hierarchy of coupled linear PDEs are obtained for the unknown functions <f) m (y). Eq. l|9} is exactly equivalent 
to the master equation and physical considerations will next be used in finding good-quality approximate solutions. We start by 
keeping only the first term in Eq. (|9jl, thus, ignoring the R — > R* dynamics, and then incorporate back the omitted terms in an 
effective way. 



A. Time-scale separation modulates noise propagation 

The first term on the right hand side of Eq. (|9} describes the birth and death of protein A and the remaining terms describe R* 
dynamics. If the R* reaction is ignored, the hierarchy of PDEs become uncoupled. We obtained an exact analytical solution for 
the resulting PDEs using the method of characteristics: 

*£?(») = </Vo[l + ( T ^— + T^-e-^A {y _ i r (10) 
\ A + fim A + fim J 

where the number of ^4's was taken to be N at t — and cj) m .o is a constant, representing the probability of having exactly m 
i?*'s. If, for example, the number m of R* is fixed at a particular value in, then 0m, o = 1 j 4>m = for to ^ fa. The obtained 
generating function indicates that the A distribution is binomial. Note that the relaxation rate is A + /ito, depending on both A 
and [i. The solution further simplifies in the limit of long times (t — > oo). 

^'-Ml + T^fe-lf , (ID 

A + [im, 

which is the generating function for the stationary distribution of A. 

In the real 2-step biochemical cascade, the number of R*'s is, of course, fluctuating. However, Eq. dlOi . with m concentrated 
at to, still constitutes a good approximation when either of the following conditions is satisfied: (i) R* is characterized by a sharp 
distribution centered at to, which is often the case when the number of R*'s is large; (ii) the reaction rates for R* birth and death 
are much larger than those for A. For a cascade that satisfies condition (i), the linear noise approximation might be applicable. 
However, our solution, based on the generating function approach, provides the full probability distribution in an analytical form. 
When condition (ii) is satisfied, the A — > A* reaction only "sees" an average number (to) of R*. In this case, our solution is 
simpler and, perhaps more convenient to use, than the ft— expansion solution. Analysis of the cascade dynamics for case (ii) 
using Eq. JlQi suggested a possible mechanism for noise filtering. Even in the case of broad or irregular distribution for R*, if 
the fluctuations around the average are fast, the distribution of A is still well approximated by Eq. dlOi . being well-peaked at the 
average for large N. 

As an example, we take the reaction rate parameter values (g,k,fj,,X) — (20,10,0.004,0.03) and initial condition 
(Nr, Nr* , Na, Na*) = (100,0,100,0). The evolution of the first two moments of Na* as computed from our approxi- 
mate solution Eq.^3 fi— expansion solution Eqs l32l3 J1 and exact numerical results are shown in Fig.0t,b. Unlike the practical 
implementation of the 51— expansion, Eq.[K)]also directly gives the time evolution of the full probability distribution for A* pro- 
teins. A time slice of the probability distribution of A* is shown in Fig.|5^. Overall, a remarkable agreement is achieved between 
the approximate analytical results and exact numerical calculations. Also shown in the figure is the distribution computed from 
the Langevin equation (dashed line), which is characterized by an average noticeably larger than the exact average. Even though 
the magnitude of R* fluctuations is the same for the cascade parameters used in Figs. 1314 l5l the A* fluctuations are dramatically 
attenuated for the case demonstrated in Fig.0J) and|5t, compared with Fig.|5] In this case the R — > R* reactions are much faster 
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FIG. 4: Time evolution of the averages ((a) and (c)) and variances ((b) and (d)) for A* obtained from three different calculations: the ap- 
proximate solution Eq. <10> (circles), the exact solution (solid line computed from Eq. Q) and the il-expansion Eq. I33t (dashed line prob- 
ably overlapped by the solid line). The initial condition is (Nr,Nr* ,Na,Na*) = (100,0, 100,0). (a) and (b) have parameter values 
(g, k, n, A) = (20, 10, 0.004, 0.03); (c) and (d) have parameter values (g, k, fi, A) = (0.02, 0.01, 0.2, 1.5). 



than the A — > A* reactions, thus, the R — > R* noise is averaged out and only internal A — > A* noise remains. If we think 
of this two-step cascade as an element of a longer signaling pathway, the relatively large fluctuations of the upstream reaction 
(R — > R*) become attenuated downstream (A — > A*). Thus, stacking of a slow downstream reaction after a fast upstream 
reaction provides a general mechanism for noise attenuation in biochemical signaling networks. 

In the opposite limit of the reaction rate of R* being much slower than that of A, we can take into account the R* dynamics 
by allowing <t> m ft in Eq. dlOl to change slowly with time according to Eq. @. Thus, we substituted 

c/> m ,o = exp(-(l - e- kt )g/k)(l - e~ kt ) m {^) m / m\ . (12) 

into Eq. dlOi . This is a valid approximation in the limit of slow R* dynamics. Using this substitution, we obtained an approximate 
solution for ^ (x, y) as an infinite sum over the R* protein number m. Since the coefficients of x m decay very fast with increasing 
in, we can safely truncate it to a finite sum. For many reaction rates in this regime, the obtained distribution is broad. 

A set of reaction rates, which corresponds to a slow upstream reaction and a fast downstream reaction (g,k,/j,,X) = 
(0.02,0.01,0.2,1.5), was used to compare our analytical calculations with exact numerical results (Fig. @];,d and Fig. |5J>). 
The evolution of both the first moment (Fig.|4];) and the variance (Fig. 01) obtained from our analytical treatment agrees well 
with the exact numerical one, being more accurate than the Sl-expansion. Furthermore, we used Eq.|8] Eq.^ll an d Ea.ll2lto 
compute the full probability distribution at t = 200 (Fig.|3J)), which is impossible to obtain analytically using the f2-expansion 
approach. All the nuances of the complicated distribution are accurately captured by our approximate analytic solution. The 
distribution computed from the Langevin equation, on the other hand, which is shown as a dashed line in Fig.|5j>, is characterized 
by a single broad peak. The white noise terms in the Langevin equation obviously smear out the peaks and, thus, are not good 
models for the underlying stochastic dynamics. At long time limits the distribution becomes uni-modal, but still wide (data not 
shown). 
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(a) (b) 

FIG. 5: Probability distributions of A* computed from the approximate solution Eg. 1101 (circles). Langevin equation (dashed line) and the 
exact solution Eq.\^(solid line) with initial condition (Nr, Nr* , Na, Na*) = (100,0, 100,0). (a) The distribution P{Na*) at t = 60 
with parameter values (g, k, /i, X) = (20, 10, 0.004, 0.03). (b) The distribution P(Na* ) at t = 200 with parameter values (g, k, fi, A) = 
(0.02,0.01,0.2,1.5). 




FIG. 6: (a) Averages computed from Gillespie simulation (solid line) and from the approximate solution Eq. |8j (circles). One typical Gillespie 
trajectory is also shown (thin solid line); (b) the variance from Gillespie simulation (solid line) and from the approximation (circles). Initial 
condition (Nr, Nr* , Na , Na* ) = (100,0,5000,0) and parameter values (g, k, fi, X) = (0.1,0.05,0.2, 1.5). 



The fluctuations in A* are much larger for a cascade where the upstream reaction is slow and the downstream reaction is 
fast (Fig. @Jl and Fig. |3J>) compared with previously considered cascades (Figs. [3] 0J> and|5^). Thus, the noise produced in 
the R — > R* reaction is retained and amplified by the next enzymatic reaction. As opposed to the previously discussed case 
of noise-attenuation, this cascade setup could be used to amplify the noise downstream. The amplification and attenuation 
of noise have been extensively discussed and experimentally tested in the linearized stochastic description of gene regulatory 
networks 55,56,58 ' 68 . Our current work emphasizes the role of discreteness in a nonlinear biochemical reaction network, using 
directly the generating function formalism to obtain analytical time-dependent probability distributions. 

Strong fluctuations occur even when the number of downstream proteins is very large. Fig. ^demonstrates a striking example. 
Although the average number of A* quickly reaches its steady-state value, which is near 1000, the fluctuation continues, such 
that a typical trajectory fluctuates with a magnitude of the same order as the average (Fig.|6^). The corresponding large variance, 
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shown in Fig.[6j), confirms this analysis. Thus, the fluctuations are much larger than expected from the usual \/N argument, 
where N is the number of proteins. Our approximate solution agrees very well with the exact solution (Fig.|6j, confirming the 
validity of our noise-amplification scheme to solve the master equation (Eqs. (|8}, dlO> and (I12» . The large A* fluctuations 
are induced by the upstream R* fluctuations, as seen from a typical numerical trajectory shown in Fig|6ji. At steady state, the 
average death/birth timescale for R* is 1/g = 10, thus, the A* trajectory exhibits a bursting behavior with a time correlation of 
about 10. 

In the photoreception cascade, the rhodopsin (Rh) activation is controlled by the incident photons which may arrive, one 
by one, within Is or even longer time interval in the single photon experiment 62 . The deactivation rate from Rh* to Rh is 
about 0.5s -1 . The activation rate of the transducing is around 120s -1 and the deactivation rate is 100s _1 . According to the 
classification developed in this paper, the initiating stage of the photoreception cascade belongs to the case of slow R* reaction 
and fast A reaction, thus, we expect an amplification of the external light stimulation^iSi^i. 

Our approximation method in this section is based on the separation of time scales for the first and the second reactions, an 
approach which is also the starting point for many solution techniques in the literature. The elimination of fast variables is among 
the most popular ones and is often used to approximately solve Fokker-Planck or Langevin equations^^. Various methods such 
as the projection operator method 32 and cumulant expansion method 31 have been developed to treat the continuous cases. To the 
lowest order approximation, the fast time scale is completely removed, while the equations for the slow variables account for 
the fast variables either through their corresponding averages 31 or through their stationary distributions^ 2 .. Here we applied the 
essence of this idea to the discrete jump process described by a master equation and obtained analytical approximations for the 
evolving probability distribution function. Similar consideration has been used to derive effective equations in gene transcription 
regulation in the large particle limit 69 , but our method applies to small particle numbers as well. In the first case considered above 
(fast A— A* reaction), we used only the average of the fast variable. In the second case (fast R—R* reaction), we first considered 
the probability distribution evolution of the fast variable and, subsequently, incorporated back the evolution of the slow variable. 
Thus, in both cases, we explicitly included both the slow and fast time scales into the final analytical expression. 



IV. SOLVING THE MASTER EQUATION WITH THE HYBRID SMOOTH PROBABILITY DISTRIBUTION METHOD 



When the probability distributions are relatively smooth, a new approximation scheme may be employed to integrate the 
generating function equation [3] To show basic idea of the method, we first implement it for the 2-step cascade discussed 
previously (Fig. Then, to demonstrate the potential for generalization, we apply this method to another enzymatic cascade, 
where the first step is a self-dimerization of receptor R instead of a simple activation. 



A. 2-step signal transduction cascade 



To treat analytically the stochastic signaling dynamics in a wider regime of parameters (for example, when upstream and 
downstream reactions have comparable rates), the neglected coupling terms in Eq|9]may be taken into account in a more sys- 
tematic way. In this section, we will reconsider the cross terms between different m's and treat them in a different way. Here, 
we take advantage of the known R* distribution from Eq.[6]and write down the following expansion for ^(x, y), 

9(x,y) = £ exp(-f (1 - e- fei ))(f ) m ?Q<l>m{y) , (13) 

ra=0 

where 4> m (y) is a time-dependent function and, according to Eq. we know that 4> m (l) = (1 — e~ kt ) m . Substituting this 
form of expansion of into Eq. l|3} and comparing the coefficient of x m on both sides of the equation, we have 



deb d d 

= (1 - y)(nm- XN + \y—)(j> m + km{4> m -i - 0m) + ff(</>m+i - <f> m ) + ge~~ kt (f> m . (14) 

ot ay oy 

This equation describes the time evolution of (f> m (y), which is coupled to the neighboring functions 4> m -i and 4> m +i . In previous 
discussions, we neglected these couplings first and only later incorporated them back in an effective way, justified under certain 
conditions. Here, we present an approach which directly takes into account these couplings in an approximate manner. When 
the probability distribution of i?*is smooth, the following approximation, 

-77"— ~ 4>m+l - 4>rn ~ 4>m ~ 4>m-l , (15) 

am 

may be used to uncouple the PDEs in Eq. MAI : 
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= (1 - tf )(/*m£ -XN + Xy^-)4> m + {g km)^ + ge^^ , (16) 
at ay ay am 

We solved the resulting PDEs by the method of characteristics. Eq. d!5l > is satisfactory when the profile of <f> m between to and 
m + 1 can be reasonably approximated by a straight line segment. This works well if either of the two conditions holds: (1) 
as mentioned above, the distribution profile is smooth so that the higher-order derivatives can be ignored, or (2) the number of 
R* is large such that an increment by one particle may be treated as small. In addition, it is possible to improve this solution by 
making a higher-order approximation to the difference in Eq. dl5l . For example, an exact expression can be obtained through 
the Kramers-Moyal expansion^ 



E 



I d< 

II dm 1 



where I runs from 1 to oo. In the standard derivation of the Fokker-Planck equation^iS, terms up to the second order (I = 2) 
are retained. In some sense, our "smooth distribution" method is a hybrid distribution function - generating function scheme, 
where in the yet to be determined functions <f> m (y), subscript m is related to the R* particle number (distribution for to), while 
y is a formal variable related to the generating function for the A* particle number. However, solving equations containing 
higher-order derivative terms quickly becomes cumbersome. Here, we only keep the first order term, which allows Eq. (II 61 to 
be solved with the following set of characteristic equations 

y = (y- l)(/J,m + \y) 
ra = mk — g 

4> m = \N{y-l)cp m +ge~ kt 4> m . (17) 

The first two equations in Eq. ( 1171 define the characteristic curves and the third equation shows how <p m changes along this 
curve. The dynamics on each curve is self-contained and independent of each other, which is the consequence of neglecting 
higher order terms in Eq. J15I . Eq. (I17> were exactly solved, resulting in 

<P rn (y) = ^>(*b)(l + ^(y - 1)) N exp[f (1 - e- kt )} , (18) 



where 



N 1 



n= Z0 

z = + \p 

y- i 

p(t) = f e'^ds 
Jo 

I(t) = (A+^+^(m -|)(e fct -l). 

a nm is the initial probability of having to R*'s and n A's and too is an intermediate variable which after all the integrations and 
differentiations are done will be replaced by to by the following substitution, 

m -| = (TO-|) e - fc *. (19) 

In general, the integration to obtain p(t) can not be carried out in a closed form, thus, approximate or numerical treatment is 
needed. However, the analytic structure of the overall solution is transparent. The last exponential factor in Eq. Jl 8b describes 
the R* reaction and the first two factors describe the A — > A* reaction. The full generating function *&(x, y) is given by Eq. <ll 31 
with 4> m (y) given by Eq. Jl 81 . Once ^(x, y) is known, all the statistical quantities are easily computable. 
As an example, we consider the following initial distribution 

'(- + 1)", (20) 
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which corresponds to starting the cascade dynamics with TV A's and a small number of R*, distributed exponentially. The 
corresponding solution is given by 



where 



P(t) 
P(t) 



<t>m{y) = (1 + ^±±(V - 1)) N exp(-me- fet ) , 
P 



exp((A+^)t+^(m-|)(l- e - fci )) 



(21) 



exp((A + -r)s + — (m - 



When x = 1, if the following approximation is used 

4> m = exp(~me- kt ) » (1 - e^')" 1 , 



(22) 



(23) 



then the R* distribution is recovered. Furthermore, Eq. d2 1 1 with Eq. (I23> substituted in conserves the total probability, i.e., 
^(1,1) = 1. Thus, we use Eq. (1231 in the calculations described below. 

To gauge the effectiveness of the approximate analytical solution, we take N — 100 and the R* distribution truncated at 
m — 30. Firstly, we evaluate Eq. d2 II with (g ,k , fx , A) = (10 , 1 , 0.01 ,0.1). A* distribution computed from Eq. Jl 31 at t — 60 
matches quite well with the exact solution as shown in Fig. Qi. The approximate distribution is a little narrower than the exact 
one due to the omission of the higher-order derivative terms in Eq. J15I . Secondly, we carried out similar calculations with 
(g ,k ,fi , A) = (0.2 , 0.1 , 0.02 , 0.2), as shown in Fig. 03. Although the distribution at t = 60 from our approximation agrees 
quite well with the exact solution for large Na*, near the left boundary there is a clear discrepancy, due to the non-smoothness 
of the distribution at the minimum particle number, since the number of A* cannot be negative. 




100 




100 



(a) The distributions at t = 6 



(b) The distributions at t = 60 



FIG. 7: The probability distributions of A* at t = 60 computed from the approximate solution Eq. 12 1 1 (circles) and the exact solution Eq. 
(solid line) with initial condition (Nr, Nr* , Na,Na* ) = (100, 0, 100, 0) and parameter values: (a) (g , k , fi , A) = (10 , 1 , 0.01 ,0.1); (b) 
(g ,k,fi,X) = (0.2 , 0.1 , 0.02 , 0.2). 



B. Receptor self-dimerization introduces additional nonlinearity 



Even if the exact solution is not available for the first reaction, we can still apply the above approximation as long as the 
distribution is smooth. Consider the dimerization reaction shown in Fig. [8] Compared to the previously discussed 2-step cascade 
(Fig. the first reaction is replaced by a self-dimerization process with rate g. This dimerization activation is quite common in 
signal transduction and gene regulatory networks2iZiS. Although it is possible to derive an analytical solution for the isolated 
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signal * * decay n n 
R + R ► R, R, — ► R + R 



A n * catalysis A * n * A » decay 
A+R, 7 - ►A+R, A ►A 



FIG. 8: An inactive receptor R, when activated by dimerization, activates downstream protein A. 



first step in the cascade, i.e. the self-dimerization process 73,74 , the expression is in a series form and will not be used in our 
approximation scheme. 

Again, let P(m, n) denote the probability of having m R^'s and n A's. The master equation is then, 

^(m,n) = | ((M - 2m + 2)(M - 2m + \)P(m- l,n) - (Af - 2m){M - 2m-l)P{m,n)) 

+k ((m + l)P(m + 1, n) — mP{m, n)) + /i ({n + l)mP(m 1 n + 1) — mnP(m, n)) 

+X ((AT - n + l)P(m, n - 1) - (JV - u)P(m, n)) , (24) 

where the first term describes the dimerization reaction. The corresponding generating function &(x, y) satisfies, 



f - f (l -l)(M(M-l) + 4^-2(2M-3,^)* 



3* a 2 * d 

+k(l - x)— + - + A(y - 1) /V - y— )* . (25) 

ox oxoy oy 

Note that a new second-order derivative, d 2 ^ / dx 2 , appears, compared with the previously considered generating function PDE 
(Eq. (|3jl). If we expand the generating function in the form of Eq. (|8jl, a series of PDEs for <p m 's are derived. Similar to the 
simple 2-step cascade, after taking the continuous limit approximation, we get 

+ 7^- (km - -(M(M + 4m(m - 1) - 2(2M - 3)mj) (f> m , (26) 
dm V 2 / 



where the total probability is conserved since 



d 
df 



dm 4> m (t, 1) = . 



Eq. J26l > are also readily solved analytically by the characteristic method. 

The approximate distributions compared with those computed from Gillespie simulation are displayed in Fig. [9] at an early 
time t = 3 and a later time t — 60. They agree quite well. The approximate and exact distributions at other times are in good 
agreement as well (data not shown). 

Overall, the presented method may be used to obtain the long time evolution of the stochastic signaling dynamics. We 
designed the hybrid scheme to treat the second order derivative terms in the generating function equation with the characteristics 
method. When the probability distribution profiles of all species are relatively smooth the method gives quantitatively accurate 
results. Although, this condition is most commonly satisfied when protein numbers are large, it is also applies to systems with 
smaller protein numbers with certain constraints on the relative reactions rates. Generalizing this method to treat larger cascades 
is also straightforward. If all the reaction nodes in the network are linear, then there is no need to expand the generating function 
equation up to second order derivative terms, and a direct application of the characteristics method would solve the generating 
function PDE. If a small number of binary reaction nodes are present, the hybrid scheme can be applied most effectively, allowing 
to obtain approximate solutions in a manner similar to the examples that we have already considered. For cascades with many 
binary reactions, a straightforward application of the method may becomes unpractical since a summation over many indices 
is required, which would be computationally expensive. However, this difficulty may be overcome by either linearizing the 
nonessential binary terms or by approximating the summations with integrals. 
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0.1 




(a) The distribution at t=3 (b) The distribution at t=60 

FIG. 9: Probability distributions of A* computed from the approximate solution Ea. \26\ (circles) and Gillespie simulation (solid line) averaged 
over 70898 realizations. Initially (Na, Na* ) = (100, 0) and Nr* approximately assumes a Gaussian distribution exp(— (m — 2) 2 )/ ^pK. We 
use N R + 2N R * = 20 and choose parameter values (g, k, fl, xf = (0.02, 0.5, 0.02, 0.15). (a) The distribution P(N A * ) at t = 3. (b) The 
distribution P(N A ' ) at i = 60. 

V. CONCLUSIONS 

Cells live in a fluctuating environment in which signals and noise keep bombarding the cell receptorsi^SiS. Noisy signals 
propagate inside the cell via microscopic chemical reaction events. The external noise interferes with the internal biochemical 
network noise originating from the underlying fundamental randomness of these chemical reactions. Cells have evolved to 
adapt to or even exploit the seemingly deleterious effect of fluctuations on signaling dynamics within a mesoscopic size object. 
Thus, it is important to develop a qualitative picture, based on mathematical modeling of stochastic chemical kinetics, of how 
signaling networks process noisy signals. In this paper we studied the stochastic signal transduction by a simple 2-step signaling 
cascade using the master equation to describe stochastic reaction events. In agreement with previous studies, we found that 
when particle numbers are large, chemical kinetic equations provide an accurate description. However, when the number of 
proteins becomes small, a large variability among individual trajectories results, necessitating the stochastic chemical kinetics 
approach. If fluctuations are small, the commonly used linear noise approximation works well, the probability distribution being 
centered around the deterministic trajectory. But when fluctuations are large, for example, at the initiating (burst) phase of a 
signal transduction cascade, the linear noise approximation breaks down and more powerful analytical treatment of the master 
equation becomes necessary. In the small protein number regime, chemical Langevin equation does not work properly as well 
since the continuous assumption breaks down and molecular discreteness sets in the dynamics. 

Without assuming that noise is small, we directly treated the master equation with a generating function approach. Although 
the resulting PDE could not be solved exactly, we found a number of perturbative schemes, that allow us to obtain approximate 
analytic solutions for the generating function which is used to obtain the time evolution of the full probability distribution for 
all proteins. Using the analytic solution, we recovered the general mechanism for attenuating or amplifying noise in a signaling 
pathway with nonlinear reaction events: if the upstream reactions are fast and the downstream reactions slow, then upstream noise 
becomes attenuated. Conversely, if the upstream reactions are slow and the downstream reactions fast, then the upstream noise 
becomes amplified. Thus, controlling various node timescales by regulating reaction rate constants would lead to enhancing 
of the signaling cascade sensitivity and reliability by suppressing uncorrelated noise while still amplifying weak signals. This 
mechanism may be used by cells to draw useful information from a noisy environment. Furthermore, under certain conditions, 
the burst phase noise may induce macroscopic system-wide disturbance in the downstream signaling network. 

The approximation based on characteristics, presented in Section llVl can be straightforwardly generalized to a longer cascade, 
with the restriction that the protein number distributions are smooth. Yet another powerful technique to solve the master equation 
for more complex cascades is based on the variational principle 76 . The analytic solutions developed in this paper serve as a 
starting point for developing high quality time-dependent basis sets for the variational approach. A good basis set should capture 
the essential part of the system dynamics so as to make the subsequent calculations simple and effective, which will be discussed 
in more detail elsewhere. For larger signaling pathways, especially when embedded in space, the commonly used numerical 
stochastic simulations will face severe computational bottleneck. In this work we have taken a different approach, based on 
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analytically solving the master equation describing stochastic chemical kinetics, to achieve, in an efficient manner, qualitative 
and quantitative insights into stochastic signaling by biochemical reaction networks. In an ongoing work we are generalizing 
the techniques presented in this work to investigate the interplay of the noise amplification and attenuation with the complex 
dynamics generated from feedback loops in larger biochemical reaction networks. 



APPENDIX A: THE fi-EXPANSION OF THE 2-STEP CASCADE MASTER EQUATION 

If the fluctuations are relatively small, van Kampen's O-expansion may be used to account for their effect on signaling 
dynamics^. In this case it is more convenient to rewrite the master equation, given by Eq.^ m an alternative form 31 , where the 
dependence of reaction rates on the cell size, CI, is explicitly emphasized: 

P = HE- 1 -l)(N- n)P + ^(E^ 1 - l)mnP + k{E+ x - l)mP + g^E' 1 - 1)P , (27) 

where £'+ 1 /(n) = /(n + 1) and E~ l f(n) = f(n — 1) are step-up and step-down operators 31 . Although CI is equal to 1 
in our units, it is still written out to be used later as an expansion parameter. If the protein numbers are large enough such 
that the deterministic evolution serves as a good starting point, then small fluctuations are well described by the linear noise 
approximation. In particular, the variables are changed to emphasize the fluctuations around the deterministic orbits: 

to = ci^t) + ci 1/2 ^ 

n = Cltp(t) + n 1/2 ij , 
P{m,n,t) = n(^TOi), 

N = NCI. (28) 

where cj)(t) traces the deterministic path for Nr* and tp(t) for Na- The new random variables are £ and 77, that describe 
fluctuations around the average path. While averages determined from the dominant paths are proportional to Cl, fluctuations 
around these averages are only proportional to CL 1 / 2 . Since <j>(t) and tp(t) can be easily found by solving the chemical kinetic 
equation, we need to obtain an evolution equation for the probability distributions of £ and 77, II(£, rj, t). Thus, we substitute 
Eq. 12 8 1 into Eq.|^J using also the following expression for any analytic function f(n): 



-!)/(») = (±^ + 5^ + ...)/(»)■ 



Ea. l27l results in 



§ - - ^ - + - •> - ^ 

(ci<t> + n^on + 9(-n- 1/2 -^ + in-^jnn + ■■■. 

We next collect terms that are of the same order in Cl. The largest Cl 1 / 2 terms give: 

d<j> on dip on w - , , <9n , , an , , an an 

which is satisfied as we choose the dynamics of <j> and ip to follow chemical kinetics 

Tt = 9 ~ k ^ 

dip 



(29) 



(// X(N - V) - \iH> ■ (30) 



At the Cl° order, we obtain 
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which is the familiar Fokker-Planck equation, derived in a systematic way. Eq.|^is the linear noise approximation to the full 
stochastic dynamics, which is valid when the path determined by Eq.[30|is stable 31 . Eq.|^is a linear PDE that has to be solved 
numerically. However, it is possible to derive a closed set of ODEs to describe time evolutions of the moments up to any order. 
For example, the averages (first moments) satisfy 



dt 
dt 



= ~k(0 , (32) 



which are just the linearized chemical kinetics equations. Note that if initial values of (ij) and (£) are taken to be zero, then they 
remain zero for all later times, consistent with the physical significance of Eq. [30]which describes the evolution of averages of 
protein numbers. This also suggests that application of Eq.|^is based on the assumption of the validity of averaged chemical 
kinetics equations. Next, we consider three second moments, satisfying 



d(y 2 ) 
dt 

dt 

d(e) 

dt 



-2(A + [i<f>)(v 2 ) - 2^(r?0 + X(N - tp) + (i</>i> 
-\(r)0 - ^(e) - »<f>(r)0 ~ HvO 

-2k{e)+W + g, (33) 



to be solved simultaneously with Eq.|30| In the current case, Eq. l3QI32l33l mav be solved analytically but the solution of PDEI31I 
is numerically cumbersome. A practical difficulty in using the Sl-expansion approach comes from the fact that Eq.|^is a (1 + 2) 
PDE, which does not seem to be a significant simplification from the master equation, Eq. In particular, similar amount of 
numerical effort is needed to obtain solutions for Eqs.^and|2] This is part of the reason we did not try to obtain the distribution 
from the fi-expansion. Therefore, we used in the main text only the moments calculated from the i7-expansion. 
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